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In  this  short  note  we  briefly  outline  a  new  and  remarkably  fast  algorithm 
for  solving  a  a  large  class  of  high  dimensional  Hamilton- Jacobi  (H-J)  initial 
value  problems  arising  in  optimal  control  and  elsewhere  [1].  This  is  done 
without  the  use  of  grids  or  numerial  approximations.  Moreover,  by  using  the 
level  set  method  [8]  we  can  rapidly  compute  projections  of  a  point  in  Mn,  n 
large  to  a  fairly  arbitrary  compact  set  [2],  The  method  seems  to  generalize 
widely  beyond  what  will  we  present  here  to  some  nonconvex  Hamiltonians, 
state  dependent  Hamiltonians,  differential  games  and  perhaps  new  linear 
programming  algorithms. 

We  begin  with  the  HJ  initial  value  problem 


t)  +  H(Vx<p{x,t))  =  0  in  Mn  x  (0,  +oo) 
(p(x,0)  =  J(x)  Vx  E  Mn. 


(1) 


We  assume  J  :  Mn  — >  M  is  convex  and  one  coercive,  i.e.,  lim||x||2_).+0O  jrO  = 
Too,  H  :  Rn  — >  M  is  convex  and  positively  one  homogeneous  (we  sometimes 
relax  all  these  assumptions). 

A  good  example  of  this  is 


H(v)  =  ||u||2. 

1 

Here  ||u||p  =  (E™=1|uj|p)p  for  p  >  1  and  (x,v)  = 

If  we  take  for  J  a  convex  Lipschitz  function  having  the  property  that, 
for  a  convex  compact  set  of  Mn 

J(x)  <0  for  any  x  E  int  U 

<  J{x)  =  0  for  any  x  E  (f2  \  int  Q) 

J(x)  >0  for  any  x  E  (Mn  \  U). 
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We  call  this  level  set  initial  data.  Then  the  set  of  points  for  which  ip(x,t)  = 
0,  t  >  0  are  exactly  those  at  a  distance  t,  from  the  boundary  of  Q.  In  fact 
given  x  ^  fl,  then  the  closest  point  .x0pt  from  x  to  dQ  is  exactly 


Vy(iM)  . 

opt  “  l|VV(x,t)||2'  W 

To  solve  (1)  we  use  the  Hopf  formula  [5] 

ip(x,t)  =(J*  +  tH)*(x )  =  —  min{  J*(v)  +  tH(v)  —  (x,  u)}, 

v€Rn 

where  the  Fenchel-Legendre  transform  /*  :  Mn  — >  R  U  (+oo)  of  the  convex 
function  /  is  defined  by 

f*(v)=  sup  {(v,x)  -  f(x)}. 
x  S  Rn 

Moreover,  for  free  we  get  that  the  minimizer  satisfies 

arg  min{J*(u)  +  tH(v )  —  (x,  v}}  =  Vx<p(x,  t).  (3) 

i>ElRn 


whenever  is  differentiable  at  x.  Let  us  note  here  that  our  algorithm 

computes  ip(x,t )  but  also  Vxtp(x,t). 

Also,  we  can  use  the  Hopf-Lax  formula  [5,  6]  to  solve  (1). 

ip(x,t)  =  min  /  J(z)  +  tH*  (— - ^  1  (4) 

z  £  Rn  ^  \  t  /  J 


for  convex  H. 

From  (4)  it  is  easy  to  show  that  if  we  have  k  different  initial  value  prob¬ 
lems  i  =  1 , . . .  k 


^-(x,t)  +  H(yxipi(x,t))  =  0,  in  Mn  x  (0,  +oo) 
0)  =  Ji{x)  Vr  6  R" 


with  the  usual  hypotheses,  then  (4)  implies,  for  any  x  £  Mn,  t  >  0 


ipi(x,t)  =  min  <  JAz)  +tH* 
z  e  Rn  ' 


x  —  z 


So 


min  LpAx.t)  =  min  <  min  <J.j(z)+tH* 


x  —  z 
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solves  the  initial  value  problem 


^  ^(x,t)  +  H(yxip(x,t))  =  0,  in  Rn 

(p(.x,0)  =  min  Ji(x)  \/x  € 

*=!,... ,fc 


x  (0,  +oo) 


(5) 


This  means  that  if  17  =  Uj=i  fcfii,  each  17,  is  compact  and  convex  and  may 
overlap,  then  we  can  easily  compute  the  set  of  all  points  at  distance  t  from 
17  which  is  exactly  the  solution  to  (5)  where  each  J;  is  a  level  set  function 
for  17, .  Moreover,  at  every  point  x  outside  of  17  for  which  there  is  one  i  such 
that  <pi(x,  t )  <  t )  for  any  i  /  i' ,  then  the  closest  point  xopt  to  x  and  17 

is  again 

Vx(fii(x,t) 

x°vt^z  W.(:M)r 

If  there  are  several  i  for  which  <pi(x,  t)  is  the  minimum  among  all  k  of  them, 
then  V 'x<p  wiU  be  “multivalued”,  i.e.  it  will  have  jumps,  but  any  of  the  xQpt 
defined  above  will  be  a  closest  point  on  17  to  x. 

We  solve  the  optimization  problem  (3)  by  using  the  split  Bregman  algo- 


3,  9]  as 

follows 

vk+1 

=  arg  min{J*(u)  —  (x,v)  +  -\\dk  —  v  —  bk  ||}, 

(6) 

2 

dk+1 

=  arg  min  j tH(d)  +  ^\\d  -  vk+1  -  bk |||| 

(7) 

bk+i 

=  bk  +  vk+1  -  dk+l . 

(8) 

Here  the  sequences  ( dk\ both  converge  to  Vxip(x,t).  Let  us  em¬ 

phasize  again  that  our  numerical  algorithm  not  only  computes  the  solution 
<p(x,t)  but  also  computes  Vxip(x,t)  when  c p(-,t )  is  differentiable. 

Both  (6)  and  (7),  up  to  change  of  variables,  can  be  reformulated  as 
finding  the  unique  minimizer  of 

argmin  1  af(w )  -| — \\w  —  zlln 
w  [  2 

which  is  the  proximal  map  of  /.  Equation  (6)  can  be  solved  if  either  J*  or 
J  have  easily  computable  proximal  maps,  which  often  occurs,  especially  if 
one  of  them  is  smooth. 

Equation  (7)  can  be  easily  solved  if  H(d)  =  ||d||2  via  the  shrink^  operator 
defined  by 


shrink^  (z,  a)  = 


max(||z||2  -a,0) 


if 

if  2  =  0 
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and  we  have 


If  H(d)  = 


argmin  <  a||u;||2  H — II  w  —  z|||l  =  shrudsMz,  a) 

w  {  2  J 

1 1  we  use  shrink i  operator  defined  as  follows  for  any  i  =  1, . . . ,  n 

Zi  —  a  if  Zi  >  a 

(shrinki(z,  a))i  =  {  0  if  \zi\  <  a 

Zi  +  a  if  Zi  <  —a 


and  we  have 

arg min  <  alkali i  H — lira  — zllil  =  shrinkUz,  a). 

w  {  2  J 

To  solve  (7)  for  more  general  H ( d )  convex  one  homogeneous  or  to  find  the 
proximal  map  for  /  of  that  type  we  use  the  fact  that  H*  is  the  characteristic 
function  of  a  closed  convex  set  C  C  Mn 


H*  =  Ic. 

By  using  the  Moreau  identity  [7]  we  realize  that  the  proximal  map  of 
H  can  be  obtained  by  projecting  onto  C.  To  do  this  projection,  we  merely 
solve  the  eikonal  equation  with  level  set  initial  data  for  C  via  split  Bregman 
as  above  in  (6),  (7),  (8)  with  H(d)  =  ||d||2-  This  is  easy  using  the  shrink2 
operator.  We  then  use  (2)  to  obtain  the  projection  and  repeat  the  entire 
iteration. 

Numerical  experiments  on  an  Intel  Laptop  Core  i5-5300U  running  at 
2.3  GHz  are  now  presented.  We  consider  diagonal  matrices  D  defined  by 
Da  =  1  +  ^  for  i  =  l,...,n.  We  also  consider  matrices  A  defined  by 
An  =  2  for  i  =  1, . . .  ,n  and  Aij  =  1  for  i,j  =  1, . . . ,  n.  Table  1  presents  the 
average  time  (in  seconds)  to  evaluate  the  solution  over  1,000,000  samples 
(x,  t)  uniformly  drawn  in  [—10, 10]n  X  [0, 10].  The  convergence  is  remarkably 
rapid:  10~6  to  10“8  seconds  on  a  standard  laptop,  per  function  evaluation. 
Figure  1  depicts  2-dimensional  slices  at  different  times  for  the  (H-J)  equation 
with  a  weighted  Hamiltonian  H  =  ||D  •  ||i,  initial  data  J  =  ^||  •  [||  and 
n  =  8. 
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n 

\\yh 

Wvh 

IMloo 

WvWd 

IlylU 

4 

6.36e-08 

1.20e-07 

2.69e-07 

7.00e-07 

8.83e-07 

8 

6.98e-08 

1.28e-07 

4.89e-07 

1.07e-06 

1.57e-06 

12 

8.72e-08 

1.56e-07 

7.09e-07 

1.59e-06 

2.23e-06 

16 

9.24e-08 

1.50e-07 

9.92e-07 

2.04e-06 

2.95e-06 
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(c)  (d) 

Figure  1:  Evaluation  of  the  solution  X2,  0,  0, 0, 0, 0, 0)t,  t)  of  the  HJ- 

PDE  with  initial  data  J  =  ^||  ■  |||  and  Hamiltonian  H  =  || D  ■  ||i  for 
(, xi,X2 )  E  [ — 20,  20]2  for  different  times  t.  Plots  for  t  =  0,3,  5,  8  and  re¬ 
spectively  depicted  in  (a),  (b),  (c)  and  (d).  The  level  lines  multiple  of  10 
are  superimposed  on  the  plots. 
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